library(rdrobust)
library(rdd)
library(tidyverse)
library(xtable)

setwd('~/Dropbox/Ideas/Women corruption/Replication_files')
#Loading main datasets for 2007 and 2011 thresholds
m07<-read.csv('m07.csv',stringsAsFactors = F)
m11<-read.csv('m11.csv',stringsAsFactors = F)

#Loading dataset with list of elected candidates in 2007 and 2011
c07<-read.csv('candidates07.csv', stringsAsFactors = F)
c11<-read.csv('candidates11.csv', stringsAsFactors = F)

# # # # ## # # # ## # # # ## # # # ## # # # ## # # # #
# Creating outcomes for difference-in-discontinuities
# # # # ## # # # ## # # # ## # # # ## # # # ## # # # #

# M07 dataset (5K threshold)  
#Revealed corruption 07-11
m07$corr7_11_diff<-m07$corr0711-m07$corr0307

#Urban dev 07-11
m07$nuc37<-ifelse(is.infinite(m07$nuc37),NA,m07$nuc37)
m07$urb7_11_diff<-m07$nuc711-m07$nuc37

#Outcomes to measure experience effects in m07 dataset (07 baseline)
m07$corr11_15_diff<-m07$corr1115-m07$corr0711
m07$UC11_15_diff<-m07$nuc1115-m07$nuc711

#Outcomes for experience effects with 2003 baseline
m07$corr11_15_diff_03baseline<-m07$corr1115-m07$corr0307
m07$UC11_15_diff_03baseline<-m07$nuc1115-m07$nuc37


# M11 dataset (3K threshold)  
#Revealed corruption 11-15
m11$corr_11_15_diff<-m11$corr1115-m11$corr0711
# With 2003 baseline
m11$corr0307<-NA
for (i in unique(m11$muni_code)){
  m11$corr0307[m11$muni_code==i]<-m07$corr0307[m07$muni_code==i]}
m11$corr_11_15_03baseline<-m11$corr1115-m11$corr0307


#Urban dev 11-15
m11$UC_11_15_diff <- m11$nuc1115 - m11$nuc711
# With 2003 baseline
m11$nuc37<-NA
for (i in unique(m11$muni_code)){
  m11$nuc37[m11$muni_code==i]<-m07$nuc37[m07$muni_code==i]}
m11$UC_11_15_03baseline<-m11$nuc1115-m11$nuc37



# # # # #
# FILTERS
# # # # #
m07$popdiff<-m07$pop-5000
m11$popdiff<-m11$pop-3000

#Municipalities where main party doesn't have voluntary party quotas
partyQ<-c("PSOE","IU","PSC","ICV","ERC","BNG","CC")
complier_m07<-m07[(m07$MainParty%in%partyQ)!=T,]
complier_m11<-m11[(m11$MainParty%in%partyQ)!=T,]
noncomplier_m07<-m07[(m07$MainParty%in%partyQ)==T,]
noncomplier_m11<-m11[(m11$MainParty%in%partyQ)==T,]
#Municipalities where share of Women candidates prior to quota adoption is already at quota levels 
compliershare_m07<-m07[m07$Wshare03<.4,]
compliershare_m11<-m11[m11$Wshare07<.4,]
noncompliershare_m07<-m07[m07$Wshare03>=.4,]
noncompliershare_m11<-m11[m11$Wshare07>=.4,]

# # # # #
# FUNCTIONS
# # # # #
sct_plots<-function(outcome = 'NA',cutoff=NA, dt=NA, label_y=NA,y_lim=NA){
  plot(dt[dt$pop<10000,outcome]~dt$pop[dt$pop<10000],pch=3,col="gray80",
       main='', xlab='Population',ylab=label_y, ylim = y_lim)
  f.l <- loess(dt[dt$pop<cutoff,outcome]~dt$pop[dt$pop<cutoff],span=.9)
  i.l <- order(dt$pop[dt$pop<cutoff])
  lines(f.l$x[i.l],f.l$fitted[i.l], lwd=2, col="red")
  f.h <- loess(dt[dt$pop>cutoff & dt$pop<10000,outcome]~dt$pop[dt$pop>cutoff & dt$pop<10000],span=.9)
  fih <- order(dt$pop[dt$pop>cutoff&dt$pop<10000])
  lines(f.h$x[fih],f.h$fitted[fih], lwd=2, col="red")
  abline(v=cutoff)
  grid(NA, NULL)
  
}

########################################
#%#%: FIGURE 1
########################################
## Proportion of (elected) female councilors in 2007 above/below 5k threshold
sct_plots(outcome = 'Wshare', cutoff = 5000, dt = m07, label_y = '% Women elected',y_lim=c(0,.8))
## Proportion of (elected) female councilors in 2011 above/below 3k threshold
sct_plots(outcome = 'Wshare', cutoff = 3000, dt = m11, label_y = '% Women elected',y_lim=c(0,.8))



########################################
#%#%: Table 2 - Immediate effects 2007 
########################################
#Revealed Corruption - Full sample
summary(rdrobust(y=m07$corr7_11_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))
#Urban Development - Full sample
summary(rdrobust(y=m07$urb7_11_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 covs=m07$nvac7[m07$pop>1000]+m07$opland07[m07$pop>1000],
                 all=TRUE))


########################################
#%#%: FIGURE 3 - Immediate effects 2007 
########################################
#Revealed Corruption
par(mar=c(5,4,2,2))
x<-rdplot(y=m07$corr7_11_diff[m07$pop<10000],
          y.lim=c(-.4,0.4),
          x.lim=c(-1000,1000),
          x=m07$popdiff[m07$pop<10000],scale=1,
          x.label="Population (threshold=5000)",
          y.label="Change in revealed corruption (07-11)",
          title="",
          col.dots='Gray', col='blue')
x$rdplot+theme(panel.grid = element_blank(),axis.text = element_text(size = 10))

#Urban Development
y<-rdplot(y=m07$urb7_11_diff[m07$pop<10000&m07$pop>1000],
          y.lim=c(-1,1),
          x.lim=c(-1000,1000),
          x=m07$popdiff[m07$pop<10000&m07$pop>1000],nbins=c(90,50),
          x.label="Population (threshold=5000)",
          y.label="Change in urban development (07-11)",
          title="",
          col.dots='Gray', col='blue')
y$rdplot+theme(panel.grid = element_blank(),axis.text = element_text(size = 10))



########################################
#%#%: FIGURE 4 - Immediate effects 2007 by treatment dosage
########################################

# LEFT PANEL - Revealed corruption
cor_hd_quotas<-rdrobust(y=complier_m07$corr7_11_diff[complier_m07$pop>1000],
                 x=complier_m07$popdiff[complier_m07$pop>1000],
                 all=TRUE)
cor_ld_quotas<-rdrobust(y=noncomplier_m07$corr7_11_diff[noncomplier_m07$pop>1000],
                        x=noncomplier_m07$popdiff[noncomplier_m07$pop>1000],
                        all=TRUE)
cor_hd_share<-rdrobust(y=compliershare_m07$corr7_11_diff[compliershare_m07$pop>1000],
                 x=compliershare_m07$popdiff[compliershare_m07$pop>1000],
                 all=TRUE)
cor_ld_share<-rdrobust(y=noncompliershare_m07$corr7_11_diff[noncompliershare_m07$pop>1000],
                 x=noncompliershare_m07$popdiff[noncompliershare_m07$pop>1000],
                 all=TRUE)

main_coefs<-c(cor_hd_quotas$Estimate[2],cor_ld_quotas$Estimate[2],cor_hd_share$Estimate[2],cor_ld_share$Estimate[2])
main_SEs<-c(cor_hd_quotas$se[3,],cor_ld_quotas$se[3,],cor_hd_share$se[3,],cor_ld_share$se[3,])
ys <- rev(c(1.1,1.3,1.9,2.1))

par(mar=c(2,7,5,1))
plot(main_coefs, ys, xlim=c(-.6,.25), ylab='',ylim=c(0.4,2.6), yaxt='n', type='n', cex=.8,
     xlab='',axes = F)
abline(v=0,lty=3)
for (i in 1:4){
  if (i%%2==1){
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='black')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='black',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='black')
  }
  else{
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='gray')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='gray',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='gray')
  }
}
axis(2,at=c(2,1.2), labels=c('Voluntary\nparty quotas', '>40% female\ncouncillors'), las=2, cex=.8,tick = F)
axis(3,at=c(-.6,-.4,-.2,0,.2))
legend('bottomleft',c('No  ','Yes'), pch=20,lty=1,col=c('black','gray'),bty='n',ncol=2,cex=.9)
mtext(side=3, "Immediate effects on\nchange in revealed corruption", padj=-1.4)


# RIGHT PANEL - Change in Urban Land

urb_hd_quotas<-rdrobust(y=complier_m07$urb7_11_diff[complier_m07$pop>1000],
                         x=complier_m07$popdiff[complier_m07$pop>1000],
                         covs=complier_m07$nvac7[complier_m07$pop>1000]+complier_m07$opland07[complier_m07$pop>1000],
                         all=TRUE)
urb_ld_quotas<-rdrobust(y=noncomplier_m07$urb7_11_diff[noncomplier_m07$pop>1000],
                 x=noncomplier_m07$popdiff[noncomplier_m07$pop>1000],
                 covs=noncomplier_m07$nvac7[noncomplier_m07$pop>1000]+noncomplier_m07$opland07[noncomplier_m07$pop>1000],
                 all=TRUE)
urb_hd_share<-rdrobust(y=compliershare_m07$urb7_11_diff[compliershare_m07$pop>1000],
                 x=compliershare_m07$popdiff[compliershare_m07$pop>1000],
                 covs=compliershare_m07$nvac7[compliershare_m07$pop>1000]+compliershare_m07$opland07[compliershare_m07$pop>1000],
                 all=TRUE)
urb_ld_share<-rdrobust(y=noncompliershare_m07$urb7_11_diff[noncompliershare_m07$pop>1000],
                 x=noncompliershare_m07$popdiff[noncompliershare_m07$pop>1000],
                 covs=noncompliershare_m07$nvac7[noncompliershare_m07$pop>1000]+noncompliershare_m07$opland07[noncompliershare_m07$pop>1000],
                 all=TRUE)


main_coefs<-c(urb_hd_quotas$Estimate[2],urb_ld_quotas$Estimate[2],urb_hd_share$Estimate[2],urb_ld_share$Estimate[2])
main_SEs<-c(urb_hd_quotas$se[3,],urb_ld_quotas$se[3,],urb_hd_share$se[3,],urb_ld_share$se[3,])
ys <- rev(c(1.1,1.3,1.9,2.1))

par(mar=c(2,7,5,1))
plot(main_coefs, ys, xlim=c(-6,2), ylab='',ylim=c(0.4,2.6), yaxt='n', type='n', cex=.8,
     xlab='',axes = F)
abline(v=0,lty=3)
for (i in 1:4){
  if (i%%2==1){
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='black')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='black',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='black')
  }
  else{
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='gray')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='gray',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='gray')
  }
}
axis(2,at=c(2,1.2), labels=c('Voluntary\nparty quotas', '>40% female\ncouncillors'), las=2, cex=.8,tick = F)
axis(3,at=seq(-6,2,by=2))
legend('bottomleft',c('No  ','Yes'), pch=20,lty=1,col=c('black','gray'),bty='n',ncol=2,cex=.9)
mtext(side=3, "Immediate effects on\nchange in urban land", padj=-1.4)


########################################
#%#%: TABLE 3 - Immediate effects 2011
########################################
#Revealed Corruption
summary(rdrobust(y=m11$corr_11_15_03baseline,
                 x=m11$popdiff,
                 all=TRUE))
#Urban development
summary(rdrobust(y=m11$UC_11_15_03baseline,
                 x=m11$popdiff,
                 covs=m11$vacland07+m11$opland07,
                 all=TRUE))

########################################
#%#%: FIGURE 5 - Immediate effects 2011
########################################
#Revealed Corruption
x<-rdplot(y=m11$corr_11_15_03baseline[m11$pop<10000&m11$pop>1000],y.lim=c(-.4,.4),x.lim=c(-1000,1000),
          x=m11$popdiff[m11$pop<10000&m11$pop>1000],scale=1,
          x.label="Population (threshold=3000)",
          y.label="Change in revealed corruption (11-15)",
          title="",
          col.dots='Gray', col='blue')
x$rdplot+theme(panel.grid = element_blank(),axis.text = element_text(size = 10))

#Urban development
y<-rdplot(y=m11$UC_11_15_03baseline[m11$pop<10000&m11$pop>1000],y.lim=c(-1,1),x.lim=c(-1000,1000),
          x=m11$popdiff[m11$pop<10000&m11$pop>1000],scale=1,
          x.label="Population (threshold=3000)",
          y.label="Change in urban land (11-15)",
          title="",
          col.dots='Gray', col='blue')
y$rdplot+theme(panel.grid = element_blank(),axis.text = element_text(size = 10))


########################################
#%#%: TABLE 4 - Experience effects 2011
########################################
#Revealed Corruption
summary(rdrobust(y=m07$corr11_15_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))
#Urban development
summary(rdrobust(y=m07$UC11_15_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))



########################################
#%#%: FIGURE 6 - Experience effects 2011
########################################
#Revealed Corruption
x<-rdplot(y=m07$corr11_15_diff[m07$pop>1000&m07$pop<10000],
          y.lim=c(-.4,0.4),
          x.lim=c(-1000,1000),
          x=m07$popdiff[m07$pop>1000&m07$pop<10000],
          x.label="Population (threshold=5000)",
          y.label="Change in revealed corruption (07-11)",scale=1,
          title="",
          col.dots='Gray', col='blue')
x$rdplot+theme(panel.grid = element_blank(),axis.text = element_text(size = 10))

#Urban development
y<-rdplot(y=m07$UC11_15_diff[m07$pop>1000&m07$pop<10000],
          y.lim=c(-1.5,1),
          x.lim=c(-1000,1000),
          x=m07$popdiff[m07$pop>1000&m07$pop<10000],scale=1.5,
          x.label="Population (threshold=5000)",
          y.label="Change in revealed corruption (07-11)",
          title="",
          col.dots='Gray', col='blue')
y$rdplot+theme(panel.grid = element_blank(),axis.text = element_text(size = 10))




# # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # #
# # # Analyses in the supplementary appendix
# # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # #

########################################
#%#%: FIGURE A1
########################################
## Change in the proportion of (elected) female councilors between 2007 and 2003 - ie Same as above with difference in discontinuities approach
m07$W07_did<-m07$Wshare-m07$Wshare03

par(mar=c(5,4,2,2))
plot(W07_did~pop,m07[m07$pop>1000&m07$pop<10000,],pch=3,col="gray80", ylim=c(-.4,.4),
     main='', xlab='Population',ylab='% Female councilors 2007 - 2003')
f.l <- loess(W07_did~pop,m07[m07$pop>1000&m07$pop<5000,],span=.9)
i.l <- order(m07$pop[m07$pop>1000&m07$pop<5000])
lines(f.l$x[i.l],f.l$fitted[i.l], lwd=2, col="red")
f.h <- loess(W07_did~pop,m07[m07$pop>5000&m07$pop<10000,],span=.9)
fih <- order(m07$pop[m07$pop>5000&m07$pop<10000])
lines(f.h$x[fih],f.h$fitted[fih], lwd=2, col="red")
abline(v=5000)
grid(NA, NULL)

## Change in the proportion of (elected) female councilors between 2011 and 2007 - ie difference in discontinuities approach
m11$W11_did<-m11$Wshare-m11$Wshare07
plot(W11_did~pop,m11[m11$pop>1000&m11$pop<10000,],pch=3,col="gray80", ylim=c(-.4,.4),
     main='', xlab='Population',ylab='% Female councilors 2011 - 2007')
f.l <- loess(W11_did~pop,m11[m11$pop>1000&m11$pop<3000,],span=.9)
i.l <- order(m11$pop[m11$pop>1000&m11$pop<3000])
lines(f.l$x[i.l],f.l$fitted[i.l], lwd=2, col="red")
f.h <- loess(W11_did~pop,m11[m11$pop>3000&m11$pop<10000,],span=.9)
fih <- order(m11$pop[m11$pop>3000&m11$pop<10000])
lines(f.h$x[fih],f.h$fitted[fih], lwd=2, col="red")
abline(v=3000)
grid(NA, NULL)

########################################
#%#%: TABLE A1
########################################
#Proportion of female COUNCILORS - 1K window
m07_1k<-m07[abs(m07$popdiff)<=1000,]
m11_1k<-m11[abs(m11$popdiff)<=1000,]

mean(na.omit(m07_1k$Wshare[m07_1k$popdiff<0]))
length(na.omit(m07_1k$Wshare[m07_1k$popdiff<0]))
mean(na.omit(m07_1k$Wshare[m07_1k$popdiff>=0]))
length(na.omit(m07_1k$Wshare[m07_1k$popdiff>=0]))
t.test(m07_1k$Wshare[m07_1k$popdiff<0],m07_1k$Wshare[m07_1k$popdiff>0])

mean(na.omit(m11_1k$Wshare[m11_1k$popdiff<0]))
length(na.omit(m11_1k$Wshare[m11_1k$popdiff<0]))
mean(na.omit(m11_1k$Wshare[m11_1k$popdiff>=0]))
length(na.omit(m11_1k$Wshare[m11_1k$popdiff>=0]))
t.test(m11_1k$Wshare[m11_1k$popdiff<0],m11_1k$Wshare[m11_1k$popdiff>0])


########################################
#%#%: TABLE A2
#########################################Proportion of female COUNCILORS - 500 window
m07_500<-m07[abs(m07$popdiff)<=500,]
m11_500<-m11[abs(m11$popdiff)<=500,]

mean(na.omit(m07_500$Wshare[m07_500$popdiff<0]))
length(na.omit(m07_500$Wshare[m07_500$popdiff<0]))
mean(na.omit(m07_500$Wshare[m07_500$popdiff>=0]))
length(na.omit(m07_500$Wshare[m07_500$popdiff>=0]))
t.test(m07_500$Wshare[m07_500$popdiff<0],m07_500$Wshare[m07_500$popdiff>0])

mean(na.omit(m11_500$Wshare[m11_500$popdiff<0]))
length(na.omit(m11_500$Wshare[m11_500$popdiff<0]))
mean(na.omit(m11_500$Wshare[m11_500$popdiff>=0]))
length(na.omit(m11_500$Wshare[m11_500$popdiff>=0]))
t.test(m11_500$Wshare[m11_500$popdiff<0],m11_500$Wshare[m11_500$popdiff>0])

########################################
#%#%: TABLE A3: Change in % of Women councilorslors (+-1000)
########################################
mean(na.omit(m07_1k$W07_did[m07_1k$popdiff<0]))
length(na.omit(m07_1k$W07_did[m07_1k$popdiff<0]))
mean(na.omit(m07_1k$W07_did[m07_1k$popdiff>=0]))
length(na.omit(m07_1k$W07_did[m07_1k$popdiff>=0]))
t.test(m07_1k$W07_did[m07_1k$popdiff<0],m07_1k$W07_did[m07_1k$popdiff>0])

mean(na.omit(m11_1k$W11_did[m11_1k$popdiff<0]))
length(na.omit(m11_1k$W11_did[m11_1k$popdiff<0]))
mean(na.omit(m11_1k$W11_did[m11_1k$popdiff>=0]))
length(na.omit(m11_1k$W11_did[m11_1k$popdiff>=0]))
t.test(m11_1k$W11_did[m11_1k$popdiff<0],m11_1k$W11_did[m11_1k$popdiff>0])

########################################
#%#%: TABLE A4: Change in % of Women councilors (+-500)
########################################
mean(na.omit(m07_500$W07_did[m07_500$popdiff<0]))
length(na.omit(m07_500$W07_did[m07_500$popdiff<0]))
mean(na.omit(m07_500$W07_did[m07_500$popdiff>=0]))
length(na.omit(m07_500$W07_did[m07_500$popdiff>=0]))
t.test(m07_500$W07_did[m07_500$popdiff<0],m07_500$W07_did[m07_500$popdiff>0])

mean(na.omit(m11_500$W11_did[m11_500$popdiff<0]))
length(na.omit(m11_500$W11_did[m11_500$popdiff<0]))
mean(na.omit(m11_500$W11_did[m11_500$popdiff>=0]))
length(na.omit(m11_500$W11_did[m11_500$popdiff>=0]))
t.test(m11_500$W11_did[m11_500$popdiff<0],m11_500$W11_did[m11_500$popdiff>0])


########################################
#%#%: TABLE A5: Change in % of Women councilors (+-1000/2000) around 3k threshold in 2007 (placebo)
########################################
#+/- 1000
m07_above3k<-m07[m07$pop>=3000 & m07$pop<4000,]
m07_below3k<-m07[m07$pop<3000 & m07$pop>2000,]
t.test(m07_below3k$Wshare,m07_above3k$Wshare)
nrow(m07_below3k)
nrow(m07_above3k)

#+/- 2000
m07_above3k<-m07[m07$pop>=3000 & m07$pop<5000,]
m07_below3k<-m07[m07$pop<3000 & m07$pop>1800,]
t.test(m07_below3k$Wshare,m07_above3k$Wshare)
nrow(m07_below3k)
nrow(m07_above3k)

# # # SECTION B

########################################
#%#%: Figure B1 - 2011 results by treatment dosage
########################################

# LEFT PANEL - Revealed corruption

cor_hd_quotas<-rdrobust(y=complier_m11$corr_11_15_03baseline,
                        x=complier_m11$popdiff,
                        all=TRUE)
cor_ld_quotas<-rdrobust(y=noncomplier_m11$corr_11_15_03baseline,
                        x=noncomplier_m11$popdiff,
                        all=TRUE)
cor_hd_share<-rdrobust(y=compliershare_m11$corr_11_15_03baseline,
                       x=compliershare_m11$popdiff,
                       all=TRUE)
cor_ld_share<-rdrobust(y=noncompliershare_m11$corr_11_15_03baseline,
                       x=noncompliershare_m11$popdiff,
                       all=TRUE)

main_coefs<-c(cor_hd_quotas$Estimate[2],cor_ld_quotas$Estimate[2],cor_hd_share$Estimate[2],cor_ld_share$Estimate[2])
main_SEs<-c(cor_hd_quotas$se[3,],cor_ld_quotas$se[3,],cor_hd_share$se[3,],cor_ld_share$se[3,])
ys <- rev(c(1.1,1.3,1.9,2.1))

par(mar=c(2,7,5,1))
plot(main_coefs, ys, xlim=c(-.2,.2), ylab='',ylim=c(0.4,2.6), yaxt='n', type='n', cex=.8,
     xlab='',axes = F)
abline(v=0,lty=3)
for (i in 1:4){
  if (i%%2==1){
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='black')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='black',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='black')
  }
  else{
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='gray')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='gray',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='gray')
  }
}
axis(2,at=c(2,1.2), labels=c('Voluntary\nparty quotas', '>40% female\ncouncillors'), las=2, cex=.8,tick = F)
axis(3,at=seq(-.2,.2,by=.1))
legend('bottomleft',c('No  ','Yes'), pch=20,lty=1,col=c('black','gray'),bty='n',ncol=2,cex=.9)
mtext(side=3, "Immediate effects on\nchange in revealed corruption", padj=-1.4)


# Right PANEL - Urban development
urb_hd_quotas<-rdrobust(y=complier_m11$UC_11_15_03baseline,
                        x=complier_m11$popdiff,
                        covs=complier_m11$vacland07+complier_m11$opland07,
                        all=TRUE)
urb_ld_quotas<-rdrobust(y=noncomplier_m11$UC_11_15_03baseline,
                        x=noncomplier_m11$popdiff,
                        covs=noncomplier_m11$vacland07+noncomplier_m11$opland07,
                        all=TRUE)
urb_hd_share<-rdrobust(y=compliershare_m11$UC_11_15_03baseline,
                       x=compliershare_m11$popdiff,
                       covs=compliershare_m11$vacland07+compliershare_m11$opland07,
                       all=TRUE)
urb_ld_share<-rdrobust(y=noncompliershare_m11$UC_11_15_03baseline,
                       x=noncompliershare_m11$popdiff,
                       covs=noncompliershare_m11$vacland07+noncompliershare_m11$opland07,
                       all=TRUE)


main_coefs<-c(urb_hd_quotas$Estimate[2],urb_ld_quotas$Estimate[2],urb_hd_share$Estimate[2],urb_ld_share$Estimate[2])
main_SEs<-c(urb_hd_quotas$se[3,],urb_ld_quotas$se[3,],urb_hd_share$se[3,],urb_ld_share$se[3,])
ys <- rev(c(1.1,1.3,1.9,2.1))

par(mar=c(2,7,5,1))
plot(main_coefs, ys, xlim=c(-3,1), ylab='',ylim=c(0.4,2.6), yaxt='n', type='n', cex=.8,
     xlab='',axes = F)
abline(v=0,lty=3)
for (i in 1:4){
  if (i%%2==1){
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='black')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='black',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='black')
  }
  else{
    segments(main_coefs[i]-1.96*main_SEs[i],ys[i], main_coefs[i]+1.96*main_SEs[i],ys[i], col='gray')
    segments(main_coefs[i]-1.65*main_SEs[i],ys[i], main_coefs[i]+1.65*main_SEs[i],ys[i], col='gray',lwd=2)
    points(main_coefs[i], ys[i], pch=20, col='gray')
  }
}
axis(2,at=c(2,1.2), labels=c('Voluntary\nparty quotas', '>40% female\ncouncillors'), las=2, cex=.8,tick = F)
axis(3,at=c(-3,-2,-1,0,.5))
legend('bottomleft',c('No  ','Yes'), pch=20,lty=1,col=c('black','gray'),bty='n',ncol=2,cex=.9)
mtext(side=3, "Immediate effects on\nchange in urban land", padj=-1.4)


# # # SECTION C

########################################
#%#%: Table C1 - Reelection rates by gender around thresholds (first time incumbents)
########################################
c07ft<-c07[c07$first_time==1,]
c11ft<-c11[c11$first_time==1,]
#Difference among women around threshold 2011
t.test(c07ft$surv11[c07ft$pop>=4000&c07ft$pop<5000 & c07ft$sexo=='M'], c07ft$surv11[c07ft$pop>=5000&c07ft$pop<6000 & c07ft$sexo=='M'])
#Difference among men around threshold 2011
t.test(c07ft$surv11[c07ft$pop>=4000&c07ft$pop<5000 & c07ft$sexo=='H'], c07ft$surv11[c07ft$pop>=5000&c07ft$pop<6000 & c07ft$sexo=='H'])
#Difference among women around threshold 2015
t.test(c11ft$surv15[c11ft$pop>=2000&c11ft$pop<3000&c11ft$sexo=='M'], c11ft$surv15[c11ft$pop>=3000&c11ft$pop<4000&c11ft$sexo=='M'])
#Difference among men around threshold 2015
t.test(c11ft$surv15[c11ft$pop>=2000&c11ft$pop<3000&c11ft$sexo=='H'], c11ft$surv15[c11ft$pop>=3000&c11ft$pop<4000&c11ft$sexo=='H'])


########################################
#%#%: Table C2 - Reelection rates by gender around thresholds (all incumbents)
########################################
#Difference among women around threshold 2011
t.test(c07$surv11[c07$pop>=4000&c07$pop<5000&c07$sexo=='M'], c07$surv11[c07$pop>=5000&c07$pop<6000&c07$sexo=='M'])
#Difference among men around threshold 2011
t.test(c07$surv11[c07$pop>=4000&c07$pop<5000&c07$sexo=='H'], c07$surv11[c07$pop>=5000&c07$pop<6000&c07$sexo=='H'])
#Difference among women around threshold 2015
t.test(c11$surv15[c11$pop>=2000&c11$pop<3000&c11$sexo=='M'], c11$surv15[c11$pop>=3000&c11$pop<4000&c11$sexo=='M'])
#Difference among men around threshold 2015
t.test(c11$surv15[c11$pop>=2000&c11$pop<3000&c11$sexo=='H'], c11$surv15[c11$pop>=3000&c11$pop<4000&c11$sexo=='H'])


# # # SECTION E

########################################
#%#%: Table E1 - Effect of 2007 quota adoption on predetermined vars
########################################

m07pred<-m07[c('corr0307','Area_hec', 'opland07', 'PP03', 'PSOE03', 'ec', 'const','tur', 'ind', 'com', 'bancos')]
m07pred$Area_hec<-m07pred$Area_hec/1000
m07pred$opland07<-m07pred$opland07/1000
m07pred$const<-m07pred$const/100

bal<-data.frame(var=NA,est=NA,std=NA,pval=NA,bw=NA)
for (i in 1:11){
  x<-rdrobust(m07pred[,i],
              m07$pop,
              c=5000)
  nm<-names(m07pred)[i]
  est<-x$coef[1]
  std<-x$se[1]
  pval<-x$pv[1]
  bw<-round(x$bws[1,1],0)
  bal<-rbind(bal,c(nm,est,std,pval,bw))
}
bal$est<-as.numeric(bal$est)
bal$std<-as.numeric(bal$std)
bal$pval<-as.numeric(bal$pval)

bal<-data.frame(bal[-1,], stringsAsFactors = F)
print(xtable(bal, digits=2),include.rownames = F)

########################################
#%#%: Table E2 - Effect of 2011 quota adoption on predetermined vars
########################################

m11pred<-m11[c('muni_code','corr0711', 'nuc711', 'Area_hec', 'opland07', 'PP07', 'PSOE07', 'ec', 'const','tur', 'ind', 'com', 'bancos')]
m11pred$Area_hec<-m11pred$Area_hec/1000
m11pred$opland07<-m11pred$opland07/1000

bal<-data.frame(var=NA,est=NA,std=NA,pval=NA,bw=NA)
for (i in 2:13){
  x<-rdrobust(m11pred[,i],
              m11$pop,
              c=3000)
  nm<-names(m11pred)[i]
  est<-x$coef[1]
  std<-x$se[1]
  pval<-x$pv[1]
  bw<-round(x$bws[1,1],0)
  bal<-rbind(bal,c(nm,est,std,pval,bw))
}
bal<-data.frame(bal[-1,], stringsAsFactors = F)
bal$est<-as.numeric(bal$est)
bal$std<-as.numeric(bal$std)
bal$pval<-as.numeric(bal$pval)

print(xtable(bal, digits=2),include.rownames = F)


########################################
#%#%: Table E3 - MAIN RESULTS WITH PRIOR TERM AS BASELINE
########################################
##### 2007 IMMEDIATE EFFECTS
#Revealed Corruption - Full sample
summary(rdrobust(y=m07$corr7_11_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))
#Urban Development - Full sample
summary(rdrobust(y=m07$urb7_11_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 covs=m07$nvac7[m07$pop>1000]+m07$opland07[m07$pop>1000],
                 all=TRUE))

##### 2011 IMMEDIATE EFFECTS WITH 2007 AS BASELINE
#Revealed Corruption
summary(rdrobust(y=m11$corr_11_15_diff,
                 x=m11$popdiff,
                 all=TRUE))
#Urban development
summary(rdrobust(y=m11$UC_11_15_diff,
                 x=m11$popdiff,
                 covs=m11$vacland07+m11$opland07,
                 all=TRUE))

##### 2011 EXPERIENCE EFFECTS WITH 2007 AS BASELINE
#Revealed corruption
summary(rdrobust(y=m07$corr11_15_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))
#Urban development
summary(rdrobust(y=m07$UC11_15_diff[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))



########################################
#%#%: Table E4 - MAIN RESULTS WITHOUT DIFF-IN-DISC
########################################
#Revealed Corruption
summary(rdrobust(y=m07$corr0711[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))
#Urban Development
summary(rdrobust(y=m07$nuc711[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 covs=m07$nvac7[m07$pop>1000]+m07$opland07[m07$pop>1000],
                 all=TRUE))

#######2011 Immediate effects (WITHOUT DIFFERENCE-IN-DISCONTINUITIES)
#Revealed Corruption
summary(rdrobust(y=m11$corr1115[m11$pop>1000],
                 x=m11$popdiff[m11$pop>1000],
                 all=TRUE))
#Urban development - Full sample
summary(rdrobust(y=m11$nuc1115[m11$pop>1000],
                 x=m11$popdiff[m11$pop>1000],
                 covs=m11$vacland07[m11$pop>1000]+m11$opland07[m11$pop>1000],
                 all=TRUE))

#######2011 Experience effects  (without difference-in-discontinuities)
#Revealed corruption
summary(rdrobust(y=m07$corr1115[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 all=TRUE))
#Urban development
summary(rdrobust(y=m07$nuc1115[m07$pop>1000],
                 x=m07$popdiff[m07$pop>1000],
                 covs=m07$nvac7[m07$pop>1000]+m07$opland07[m07$pop>1000],
                 all=TRUE))


########################################
#%#%: Figure E1 - McCrary sorting test
########################################
###############McCrary Sorting Test (Fig E1)
dev.off()
DCdensity(runvar=m07$popdiff[m07$pop<10000],cutpoint=0,bin=200,bw=1000, verbose=T)
abline(v=0)
mtext("2007 threshold (5,000)", side=1, 2.5)
DCdensity(runvar=m11$popdiff[m11$pop<10000],cutpoint=0,bin=200,bw=1000, verbose=T)
abline(v=0)
mtext("2011 threshold (3,000)", side=1, 2.5)

########################################
#%#%: Figure E2 - Bandwidth sensitivity
########################################
#07 SHORT TERM REVEALED CORRUPTION
rg_rc7<-seq(from=100, to=2500, by=100)
sens_rc7<-rep(NA,5)
for (i in rg_rc7){
  mod<-rdrobust(m07$corr7_11_diff[m07$pop>1000],
                m07$pop[m07$pop>1000],
                c=5000,
                h=i)
  ci1<-mod$ci[1,]
  ci2<-mod$ci[3,]
  est<-mod$coef[1]
  sens_rc7<-rbind(sens_rc7,c(est, ci1, ci2))  
}
sens_rc7<-sens_rc7[-1,]

# PLOT 1
par(mar=c(4,4,1,2))
plot(rg_rc7,sens_rc7[,1],type='l',
     xlab="Bandwidth", ylab='Estimate of change in revealed corruption',
     xlim=c(0,2500), ylim=c(-.8,.4),
     main="")
lines(rg_rc7,sens_rc7[,2], lty=2)
lines(rg_rc7,sens_rc7[,3], lty=2)
abline(0,0, col="Gray")
abline(v=1052,col='blue')


#07 SHORT TERM URBAN DEV
rg_uc7<-seq(from=100, to=2500, by=100)
sens_uc7<-rep(NA,5)
for (i in rg_uc7){
  mod<-rdrobust(y=m07$urb7_11_diff[m07$pop>1000],
                x=m07$popdiff[m07$pop>1000],
                covs=m07$nvac7[m07$pop>1000]+m07$opland07[m07$pop>1000],
                h=i)
  ci1<-mod$ci[1,]
  ci2<-mod$ci[3,]
  est<-mod$coef[1]
  sens_uc7<-rbind(sens_uc7,c(est, ci1, ci2))  
}
sens_uc7<-sens_uc7[-1,]

#PLOT 2
plot(rg_uc7,sens_uc7[,1],type='l',
     xlab="Bandwidth", ylab='Estimate of change in urban development',
     xlim=c(0,2500), ylim=c(-8,5),
     main="")
lines(rg_uc7,sens_uc7[,2], lty=2)
lines(rg_uc7,sens_uc7[,3], lty=2)
lines(rg_uc7,sens_uc7[,1])
abline(0,0, col="Gray")
abline(v=1113,col='blue')

#11 SHORT TERM REVEALED CORRUPTION
rg_rc11<-seq(from=100, to=2500, by=100)
sens_rc11<-rep(NA,5)
for (i in rg_rc11){
  mod<-rdrobust(m11$corr_11_15_03baseline, m11$pop,
                c=3000,
                h=i)
  ci1<-mod$ci[1,]
  ci2<-mod$ci[3,]
  est<-mod$coef[1]
  sens_rc11<-rbind(sens_rc11,c(est, ci1, ci2))  
}
sens_rc11<-sens_rc11[-1,]

# Plot 3
plot(rg_rc11,sens_rc11[,1],type='l',
     xlab="Bandwidth", ylab='Estimate of change in revealed corruption',
     xlim=c(0,2500), ylim=c(-.8,.4),
     main="")
lines(rg_rc11,sens_rc11[,2], lty=2)
lines(rg_rc11,sens_rc11[,3], lty=2)
abline(0,0, col="Gray")
abline(v=1136,col='blue')

#11 SHORT TERM URBAN DEVELOPMENT
rg_uc11<-seq(from=100, to=2500, by=100)
sens_uc11<-rep(NA,5)
for (i in rg_uc11){
  mod<-rdrobust(m11$UC_11_15_03baseline, 
                m11$pop,
                c=3000,
                h=i)
  ci1<-mod$ci[1,]
  ci2<-mod$ci[3,]
  est<-mod$coef[1]
  sens_uc11<-rbind(sens_uc11,c(est, ci1, ci2))  
}
sens_uc11<-sens_uc11[-1,]

# Plot 4
plot(rg_uc11,sens_uc11[,1],type='l',
     xlab="Bandwidth", ylab='Estimate of change in urban development',
     xlim=c(0,2500), ylim=c(-8,5),
     main="")
abline(0,0, col="Gray")
lines(rg_uc11,sens_uc11[,2], lty=2)
lines(rg_uc11,sens_uc11[,3], lty=2)
abline(v=1055,col='blue')


#EXPERIENCE REVEALED CORRUPTION
exp_rc7<-seq(from=100, to=2500, by=100)
sens_rc7<-rep(NA,5)
for (i in exp_rc7){
  mod<-rdrobust(m07$corr11_15_diff,
                m07$pop,
                c=5000,
                h=i)
  ci1<-mod$ci[1,]
  ci2<-mod$ci[3,]
  est<-mod$coef[1]
  sens_rc7<-rbind(sens_rc7,c(est, ci1, ci2))  
}
sens_rc7<-sens_rc7[-1,]

# PLOT 5
plot(exp_rc7,sens_rc7[,1],type='l',
     xlab="Bandwidth", ylab='Estimate of change in revealed corruption',
     xlim=c(0,2500), ylim=c(-.8,.4),
     main="")
lines(exp_rc7,sens_rc7[,2], lty=2)
lines(exp_rc7,sens_rc7[,3], lty=2)
abline(0,0, col="Gray")
abline(v=1740,col='blue')


#EXPERIENCE URBAN DEV
exp_uc7<-seq(from=100, to=2500, by=100)
sens_uc7<-rep(NA,5)
for (i in exp_uc7){
  mod<-rdrobust(m07$UC11_15_diff, m07$pop,
                c=5000,
                h=i)
  ci1<-mod$ci[1,]
  ci2<-mod$ci[3,]
  est<-mod$coef[1]
  sens_uc7<-rbind(sens_uc7,c(est, ci1, ci2))  
}
sens_uc7<-sens_uc7[-1,]

#PLOT 6
plot(exp_uc7,sens_uc7[,1],type='l',
     xlab="Bandwidth", ylab='Estimate of change in urban development',
     xlim=c(0,2500), ylim=c(-8,5),
     main="")
lines(exp_uc7,sens_uc7[,2], lty=2)
lines(exp_uc7,sens_uc7[,3], lty=2)
lines(exp_uc7,sens_uc7[,1])
abline(0,0, col="Gray")
lines(exp_uc7,sens_uc7[,1])
abline(v=2295,col='blue')


